{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from IPython.display import IFrame\n",
    "import pylab as plt\n",
    "%matplotlib inline\n",
    "import urllib\n",
    "# You only need this line if you haven't cloned the repo...if you have cloned, you'll already have the data\n",
    "urllib.urlretrieve('https://raw.githubusercontent.com/sdrogers/fcmlcode/master/notebooks/data/olympic100m.txt', 'olympic100m.txt')\n",
    "import numpy as np\n",
    "# If you have cloned, make sure this is pointing to the correct file, maybe ../data/olympic100m.txt ?\n",
    "data = np.loadtxt('olympic100m.txt',delimiter=',')\n",
    "x = data[:,0][:,None]\n",
    "t = data[:,1][:,None]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear regression in vector and matrix format\n",
    "\n",
    "In the previous notebook we derived the values of $w_0$ and $w_1$ that minimised the loss. If we modified our model to have a quadratic term:\n",
    "\n",
    "$$ t_n = w_0 + w_1x_n + w_2x_n $$\n",
    "\n",
    "resulting in more exciting bendy lines, we could do the same process:\n",
    "\n",
    " - Multiply out the squared loss\n",
    " - Differentiate with respect to $w_0$, $w_1$ and $w_2$\n",
    " - Set to zero\n",
    " - Solve the resulting 3 simultaneous equations\n",
    " \n",
    "Life is too short for this. So, we re-pose the problem using *vectors* and *matrices*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Definitions\n",
    "\n",
    "We'll start with the linear model again, and define the following:\n",
    "\n",
    "We first define $\\mathbf{w},\\mathbf{x}_n$ as:\n",
    "\n",
    "$$ \\mathbf{w} = \\left[\\begin{array}{c} w_0\\\\w_1\\end{array}\\right],~~\\mathbf{x}_n = \\left[\\begin{array}{c} 1 \\\\ x_n \\end{array}\\right] $$\n",
    "\n",
    "and:\n",
    "\n",
    "$$ \\mathbf{t} = \\left[\\begin{array}{c} t_1\\\\ t_2\\\\ \\vdots\\\\ t_N \\end{array}\\right] $$\n",
    "\n",
    "and:\n",
    "\n",
    "$$ \\mathbf{X} = \\left[\\begin{array}{c} \\mathbf{x}_1^T \\\\ \\mathbf{x}_2^T \\\\ \\vdots \\\\ \\mathbf{x}_N^T \\end{array}\\right] = \\left[ \\begin{array}{cc} 1 & x_1 \\\\ 1 & x_2 \\\\ \\vdots \\\\ 1 & x_N \\end{array}\\right] $$\n",
    "\n",
    "This might all seem a bit odd (particularly the $\\mathbf{X}$), but as you see in the lecture, it means we can write:\n",
    "\n",
    "$$ \\mathbf{t} = \\mathbf{X}\\mathbf{w} $$\n",
    "\n",
    "and \n",
    "\n",
    "$$ L = \\frac{1}{N}(\\mathbf{t} - \\mathbf{X}\\mathbf{w})^T (\\mathbf{t} - \\mathbf{X}\\mathbf{w}) $$\n",
    "\n",
    "which can be differentiated with respect to $\\mathbf{w}$ (as shown in lectures and sketched below) to give the following optimal value of $\\mathbf{w}$:\n",
    "\n",
    "$$ \\mathbf{w} = \\left(\\mathbf{X}^T\\mathbf{X}\\right)^{-1}\\mathbf{X}^T\\mathbf{t} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"600\"\n",
       "            height=\"300\"\n",
       "            src=\"vec_diff.pdf\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x104989210>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "IFrame('vec_diff.pdf',width=600,height=300)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now construct the various objects we need..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "X = np.hstack((np.ones_like(x),x))\n",
    "t = t # This is already a vector!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And compute $\\mathbf{w}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  3.64164559e+01]\n",
      " [ -1.33308857e-02]]\n"
     ]
    }
   ],
   "source": [
    "XX = np.dot(X.T,X)\n",
    "invXX = np.linalg.inv(XX)\n",
    "Xt = np.dot(X.T,t)\n",
    "w = np.dot(invXX,Xt)\n",
    "print w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can print the results. To do so, we need to create a matrix of test points in the same format as $\\mathbf{X}$ above (people often get stuck here!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "testx = np.linspace(1896,2012,100)[:,None]\n",
    "testX = np.hstack((np.ones_like(testx),testx))\n",
    "testt = np.dot(testX,w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x106315790>]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAAFkCAYAAABmeZIKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3Xl4VNX9x/H3dyCCQERbiygiCeIScMHglgquCNgKomg1\n1brvxiCgdQELKO4gguK+URf8qahAFSjuihFrIigSVJaoRcUdAmqN5vz+OBMZwgSSydzZ8nk9Tx4y\ndz1zSO58cu4595hzDhEREZEghJJdABEREclcChoiIiISGAUNERERCYyChoiIiARGQUNEREQCo6Ah\nIiIigVHQEBERkcAoaIiIiEhgFDREREQkMAoaIiIiEpgGBw0z62Vm081shZlVm9mAiHXNzewGM3vX\nzNaEt5lsZtvW47gXmNlyM/vRzN40s30aWjYRERFJLbG0aLQG5gMXALUnSmkFdAdGA3sBRwO7ANM2\ndkAzOx4YB4wM77cAmG1mW8dQPhEREUkR1phJ1cysGhjonJu+kW32BuYBnZxz/61jmzeBec65weHX\nBnwKTHTO3RhzAUVERCSpEtFHY0t8y8f30VaaWRbQA3ihZpnz6ed5oCAB5RMREZGANA/y4GbWArge\neNQ5t6aOzbYGmgEray1fib/tEu24vwf6AhXAT3EprIiISNPQEsgBZjvnvgn6ZIEFDTNrDjyBb804\nP5ZDsGEfkBp9gUdiLJqIiIjAicCjQZ8kkKARETI6AodupDUD4GvgV2CbWsvbsWErR40KgIcffpi8\nvLzGFTbJzj3ySO74/HMsyjoHnLftttz5r39t9BhDhgxh/PjxgZQvnage1lFdeKqHdVQXnuoBysvL\nOemkkyD8WRq0uAeNiJDRGTjEOffdxrZ3zlWZWSlwGDA9fAwLv55Yx24/AeTl5ZGfnx+voifFEcce\ny1eTJtGvunqDdTNDIf503HGbfI9t27ZN+3qIB9XDOqoLT/WwjurCUz2sJyFdD2J5jkZrM9vTzLqH\nF3UOv+5oZs2AqUA+cBKQZWbbhL+yIo7xgplF3k65GTjbzE42s12BO/FDZR+M8X2ljYuvuYab8/KY\nGQr9dp/I4UPG+Lw8ho0Zk8ziiYiINEosLRp7Ay/hPw8d/vkXAJPxz8/oH14+P7y8pq/FIcCr4WW5\n+E6gADjnHg8/M+Mq/C2U+UBf59xXMZQvrWRnZzO1pIRxI0Zw8/TptKqq4oesLA4YMICpY8aQnZ2d\n7CKKiIjErMFBwzn3ChtvCdlkK4lzrnOUZbcDtze0PJkgOzubURMmwIQJOOfwd45ERETSn+Y6STGx\nhIzCwsIASpJ+VA/rqC481cM6qgtP9ZB4jXoyaLKYWT5QWlpaqk49IiIiDVBWVkaPHj0AejjnyoI+\nn1o0REREJDAKGiIiIhIYBQ0REREJjIKGiIiIBEZBQ0RERAKjoCEiIiKBUdAQERGRwChoiIiISGAU\nNERERCQwChoiIiISGAUNERERCYyChoiIiARGQUNEREQCo6AhIiIigVHQEBERkcAoaIiIiEhgFDRE\nREQkMAoaIiIiEhgFDREREQmMgoaIiIgERkFDREREAqOgISIiIoFR0BAREZHAKGiIiIhIYBQ0RERE\nJDAKGiIiIhIYBQ0REREJjIKGiIiIBEZBQ0RERAKjoCEiIiKBUdAQERGRwChoiIiISGAUNERERCQw\nChoiIiISGAUNERERCYyChoiIiARGQUNEREQCo6AhIiIigVHQEBERkcAoaIiIiEhgFDREREQkMAoa\nIiIiEhgFDREREQmMgoaIiIgERkFDREREAqOgISIiIoFR0BAREZHAKGiIiIhIYBQ0REREJDANDhpm\n1svMppvZCjOrNrMBtdYfbWazzOyr8Po96nHMU8Lb/hr+t9rMfmho2URERCS1xNKi0RqYD1wAuDrW\nvw5cWsf6uqwC2kd8dYqhbCIiIpJCmjd0B+fcLGAWgJlZlPUPh9d1AjZYv/FDu68aWh4RERFJXanU\nR6ONmVWY2Sdm9oyZdU12gURERKRxUiVofACcDgwATsSX6w0z65DUUomIiEijNPjWSRCcc28Cb9a8\nNrMSoBw4GxhZ135Dhgyhbdu26y0rLCyksLAwoJKKiIikjylTpjBlypT1lq1atSqhZTDnGtJfs9bO\nZtXAQOfc9CjrOgHLge7OuXdjOPbjQJVz7sQo6/KB0tLSUvLz82MouYiISNNUVlZGjx49AHo458qC\nPl/Qt05iSjFmFgJ2Az6Pb3FEREQkkRp868TMWgNdWDeipLOZ7Ql865z71My2AnYAOoS32TU8OuUL\n59zK8DEmAyucc1eEX1+Jv3WyBNgS+Dt+eOu9jXlzIiIiklyx9NHYG3gJ31rhgHHh5ZNZ16HzgYj1\nNTeHRgNXhb/vCPwaccytgLvxz8/4DigFCpxzi2Mon4iIiKSIWJ6j8QobueXinJuMDx0bO8ahtV4P\nBYY2tCwiIiKS2lJleKuIiIhkIAUNERERCYyChoiIiARGQUNEREQCo6AhIiIigVHQEBERkcAoaIiI\niEhgFDREREQkMAoaIiIiEhgFDREREQmMgoaIiIgERkFDREREAqOgISIiIoFR0BAREZHAKGiIiIhI\nYBQ0REREJDAKGiIiIhIYBQ0REREJjIKGiIiIBEZBQ0RERAKjoCEiIiKBUdAQERGRwChoiIiISGAU\nNERERCQwChoiIiISGAUNERERCYyChoiIiARGQUM24JxLdhFERCRDKGgIAJWVlYwsLqZ3bi4DO3ak\nd24uI4uLqaysTHbRREQkjTVPdgEk+SorKxlUUMDQ8nJGVVdjgANmT5rEoBdfZGpJCdnZ2ckupoiI\npCG1aAhjhw9naHk5/cIhA8CAftXVDCkvZ9yIEcksnoiIpDEFDWHujBn0ra6Ouq5fdTVzp09PcIlE\nRCRTKGg0cc45WldV/daSUZsBraqq1EFURERioqDRxJkZa7OyqCtGOGBtVhZmdUURERGRuiloCAf0\n78/sUPQfhVmhED0HDEhwiUREJFMoaAgXX3MNN+flMTMU+q1lwwEzQyHG5+UxbMyYZBZPRETSmIKG\nkJ2dzdSSEuYVFdEnJ4ejOnSgT04O84qKNLRVREQaRc/REMCHjVETJsCECTjn1CdDRETiQi0aKWL1\nath/f5g8GX75JbllUcgQEZF4UdBIEatXQ/v2cOqpsOuu8OCDyQ8cIiIijaWgkSK23x6eeQbKymD3\n3eG003zgeOABqKpKdulERERio6CRYvbaC55+Gt55B/bYA04/3QeO++9X4BARkfSjoJGiuneHp56C\n+fP992ecAbvsAvfdp8AhIiLpQ0Ejxe25J0ydCgsWQH4+nHkm7Lwz3Hsv/PxzsksnIiKycQoaaWKP\nPeDJJ+Hdd2GffeDss33guPtuBQ4REUldChppZvfd4fHHfeDYf38491zYaSe46y4FDhERST0KGmlq\nt93gscfgvfegoADOOw+6dIE774T//S/ZpRMREfEUNNJct24+cCxcCD17wvnn+xYOBQ4REUkFChoZ\nomtXePRReP996NXLB44uXeD22xU4REQkeRQ0MkxeHjzyCCxaBAcdBBdeCDvuCLfdBj/9lOzSiYhI\nU6OgkaF23RUeftgHjkMOgcGDFThERCTxGhw0zKyXmU03sxVmVm1mA2qtP9rMZpnZV+H1e9TzuMeZ\nWbmZ/WhmC8zsiIaWralyztW5bpdd4KGHoLwcevdeFzgmTnQKHCIiErhYWjRaA/OBC4Bon3CtgdeB\nS+tYvwEzKwAeBe4BugPPAM+YWdcYytckVFZWMrK4mN65uQzs2JHeubmMLC6msrIy6vY77wy33VbJ\neSdeTfWqqQweXE3bNivpd+BUvvwy+j4iIiKN1byhOzjnZgGzACzKfOLOuYfD6zoB9Z1vfDAw0zl3\nc/j1SDPrAxQB5ze0jJmusrKSQQUFDC0vZ1R1NYZPdLMnTWLQiy8ytaSE7OzsOve5tbqapezINb8O\n55+v/Y2OHb5j9NVZDB7cks03T8pbEhGRDJUqfTQKgOdrLZsdXi61jB0+nKHl5fQLhwzwia5fdTVD\nyssZN2LEJvfpwlIe4HQ+YFcO/OVZrrgii9xcGD8efvghke9GREQyWaoEjfbAylrLVoaXSy1zZ8yg\nb3V11HX9qquZO316vffpwlL+zWkUbHcof/4zXHIJdO4M48YpcIiISOM1+NZJAtXcEajTkCFDaNu2\n7XrLCgsLKSwsDLJcSeWco3VVVZ33pAxoVVWFc46aO1v12WdrlnLvvY7hw41rr4XLLoMbb/TB47zz\noHXrAN6MiIgEasqUKUyZMmW9ZatWrUpoGVIlaHwBbFNrWTs2bOVYz/jx48nPzw+sUKnIzFiblYUj\negcYB6zNyiKy+0xD9unc2c8MO3w4XHstXH75usBx/vkKHCIi6STaH99lZWX06NEjYWUI+tZJvUad\nACXAYbWWHR5eLrUc0L8/s0PR/+tmhUL0HDBgg+UN3Sc3F+65Bz76CI4+Gq64wi+78UZYs6bx70FE\nRJqGWJ6j0drM9jSz7uFFncOvO4bXb2VmewLd8H9A7xpev03EMSab2bURh50AHGFmQ81sFzMbBfQA\nbovxfWW0i6+5hpvz8pgZCv2W5BwwMxRifF4ew8aMics+ADk5fmbYJUvgmGNgxAgfOG64QYFDREQ2\nLZYWjb2Bd4BS/GfVOKAMGB1ePyC8fkZ4/ZTw+nMijtGRiI6ezrkSoBA4G/+MjmOAo5xzi2IoX8bL\nzs5makkJ84qK6JOTw1EdOtAnJ4d5RUVRh7bGuk+kTp38RG1LlsCxx8KVV/oQcv31UMejO0RERLCN\nPVUyVZlZPlBaWlra5PpoRBPZ8TPIfSJ98okPGffdB9nZMGwYFBX570VEJHVF9NHo4ZwrC/p8qTK8\nVRohlsDQmJABsMMOfmbYJUvg+ONh1CjfwnHttbB6daMOLSIiGURBQxqlY0eYNAmWLoXCQhg92geO\nMWMUOEREREFD4mT77f3MsEuXwokn+qCRkwNXXw0JHrItIiIpREFD4mr77eHWW33gOOkkuOYaHziu\nugq+/z7ZpRMRkURT0JBAdOgAEyfCsmVwyilw3XU+cIwercAhItKUKGhIoLbbDm65xQeO007zI1Vy\ncnznUQUOEZHMp6AhCbHttn5m2OXL4fTT/RNGO3WCf/wDvvsu2aUTEZGgKGhIQrVvDzff7Fs4zjwT\nxo71LRxXXgnffpvs0omISLwpaEhStG/vp6JfvhzOOst/n5PjH3H+zTfJLp2IiMSLgoYk1Tbb+FaN\nigo491x/eyUnx88eq8AhIpL+FDQkJbRr5/ttVFT46egnTPCB44or4Ouvk106ERGJlYKGpJQ//MHP\nDFtR4edOmTjRB47LLoOvvkp26UREpKEUNCQlbb21f/ZGRQVceKF/zHluLlx6qQKHiEg6UdCQlBYZ\nOAYP9hO55eTA3/8OX36Z7NKJiMimKGhIvTnnErpfpN//3j/OvKICLroI7rzTt3BcfDGsXBm/89RX\nIs8lIpLOFDRkoyorKxlZXEzv3FwGduxI79xcRhYXU1lZGch+mxIZOIYOhbvvdmy//c90ansv/bbL\nj9t5ognqPYmIZDTnXNp9AfmAKy0tdRKc1atXu8O7dXMzQyFXDc6BqwY3MxRyh3fr5lavXh3X/WIp\n38G7FLhCrnJb8L3bnLXuIsa5R2zbuJ6n5lyJeE8iIkErLS11gAPyXQI+s9WiIXUaO3w4Q8vL6Vdd\njYWXGdCvupoh5eWMGzEirvvFUr5LP5rHo/yDCnK4hJu4nzM4wy2l9ftnMmrIjXE5T825EvGeREQy\njYKG1GnujBn0ra6Ouq5fdTVzp0+P636NKd9WfM9oRlFBDpdxPS9xKrfcP5yLLoLPP4/vuWqL53sS\nEck0ChoSlXOO1lVVv/31XpsBraqqNugUGet+8SrfVnzPSK6ighx2anMbkyc7cnP9iJXPPovvuWrE\n6z2JiGQiBQ2JysxYm5VFXR+dDliblYXZ+h+/se4X7/K1ZRXb/34SFRXGiBHwz39C585QXAwrVsT3\nXPF6TyIimUhBQ+p0QP/+zA5F/xGZFQrRc8CAuO4XRPnatvUTtVVU+H8ffhh23NE/dfS//43vuURE\nJIpE9DiN9xcadZIQNSMtnqs10uK5eo46aeh+iSjfqlXOjRnj3O9+59xmmzl3/vnOffJJMOcSEUlF\nGnUiKSM7O5upJSXMKyqiT04OR3XoQJ+cHOYVFTG1pITs7Oy47peI8m2xhZ8ZdvlyGDkSHnsMunTx\nE7l9+mn860JEpKkzl4Yd2MwsHygtLS0lPz8/2cVpMpxzMfVDiHW/RJynstLPozJ2LKxeDWecAZdf\nDjvsEP9ziYikgrKyMnr06AHQwzlXFvT51KIh9RbrB2uiPpBjOU92tp8ZdvlyuPpqePJJ38Jx7rnw\n8cfxPZeISFOkoCGCDxyXXuoDx5gxMHUq7LQTnHPOxgOHiIhsnIKGSIQ2bfzMsMuX+zlVnn7at3Cc\nfbYfuSIiIg2joCESRZs2cMklPnBcdx1Mm+ZbOM48E5YtS3bpRETSh4KGyEa0bu2nol+2DK6/HmbM\ngJ139p1GFThERDZNQUOkHlq3hmHDfAvHTTfBs8/6wHHaabBkSbJLJyKSuhQ0RBqgVSsYMsS3Zowd\nC7Nmwa67wqmnKnCIiESjoCESg1at4KKLfOAYNw7+/W8fOE45BT76KNmlExFJHQoaIo2w+eZ+Ztil\nS+Hmm2HOHB84Tj4ZPvww2aUTEUk+BQ2RONh8cz8z7LJlcMst8MILkJcHf/sbfPBBsksnIpI8Choi\ncdSyJVx4oW/hmDgRXnoJunaFk06CxYuTXToRkcRT0BAJQMuWcMEFvoPorbfCK6/4wHHiiQocItK0\nKGiIBKhlSz8z7JIlfvK2117zgaOwEBYtSnbpRESCp6AhkgAtWsB55/kRKbffDnPnwm67wQknwPvv\nJ7t0IiLBUdCQJs85l7BztWjhZ4b96CO44w4oKYHdd4fjj4eFCxNWDBGRhFHQkCapsrKSkcXF9M7N\nZWDHjvTOzWVkcTGVlZUJOX+LFn5m2I8+grvugnnzfOD4y18UOEQksyhoSJNTWVnJoIICCiZNYk5F\nBdNWrGBORQUFkyYxqKAgYWEDYLPN4Kyz/DM37rkH/vMfHziOOw7eey9hxRARCYyChjQ5Y4cPZ2h5\nOf2qq7HwMgP6VVczpLyccSNGJLxMm23mZ4b98EO4914oLYU99oBBg2DBgoQXR0QkbhQ0pMmZO2MG\nfauro67rV13N3OnTE1yidbKy/MywH3wA990H8+dD9+5wzDH+exGRdKOgIU2Kc47WVVW/tWTUZkCr\nqqqEdhCNJisLTj/dP3PjgQfg3Xdhr73g6KPhnXeSWjQRkQZR0JAmxcxYm5VFXTHCAWuzsjCrK4ok\nVlaWnxl28WJ48EHfbyM/HwYOVOAQkfSgoCFNzgH9+zM7FP1Hf1YoRM8BAxJcok1r3tzPDLt4MUye\n7J+9kZ8PRx0FZWXJLp2ISN0UNKTJufiaa7g5L4+ZodBvLRsOmBkKMT4vj2FjxiSzeBvVvLmfGba8\n3AeO8nLo0QMGDPAdSEVEUo2ChjQ52dnZTC0pYV5REX1ycjiqQwf65OQwr6iIqSUlZGdnJ7uIm1QT\nOBYtgoce8p1H994b+veHt99OdulERNaxZHd6i4WZ5QOlpaWl5OfnJ7s4kuaccynTJyNWv/4Kjz0G\nV1/tQ8ef/gQjR8K++ya7ZCKSasrKyujRowdAD+dc4Ddf1aIhTV66hwyAZs38zLDvvw+PPgrLlsF+\n+/nAMW9esksnIk1Zg4OGmfUys+lmtsLMqs1sg55zZnaVmX1mZj+Y2Rwz67KJY44MHyvyS3NbijRQ\ns2Z+ZtiFC2HKFKiogP33hyOOgDffTHbpGicdW19FJLYWjdbAfOAC2HCUoJldChQB5wD7AmuB2Wa2\n2SaOuxDYBmgf/uoZQ9lEBB84TjjBD4d97DH45BMoKIC+ff1Ebuki2XPSiEjjNThoOOdmOef+4Zx7\nBqI+92gwcLVzboZzbiFwMrAdMHATh/7FOfeVc+7L8Ne3DS2biKyvWTM/M2xN4Pjvf+GPf4Q+feCN\nN5Jduo1LpTlpRCR2ce2jYWa5+NaIF2qWOedWA/OAgk3svlP4dsxSM3vYzDrGs2wiTVkotC5wPP44\nfPYZHHAAHH44vP56sksXXSrOSSMiDRfvzqDt8bdTVtZavjK8ri5vAqcCfYFzgVzgVTNrHefyiTRp\noZCfGfbdd+GJJ2DlSujVC3r3htdeS3bp1pfKc9KISP01T9B5jCj9OWo452ZHvFxoZm8BHwN/AR6o\na78hQ4bQtm3b9ZYVFhZSWFjYuNKKZLhQCI491k/W9swzcNVVcOCBcOihfljsgQcmt3wNmZMmE0YN\niQRlypQpTJkyZb1lq1atSmgZ4h00vsBfA7Zh/VaNdkC9Z2Zwzq0ysw+BjY5WGT9+vJ6jIdIIoZAP\nGwMHwvTpMGoUHHQQHHKIDxwHHZScckXOSRMtRqTanDQiqSraH98Rz9FIiLjeOnHOLceHjcNqlpnZ\nFsB+QL27nplZG2BH4PN4lk9EoguFfNgoK4Onn4bvvoODD/aB4+WXk1OmdJyTRkQ2FMtzNFqb2Z5m\n1j28qHP4dU3nzVuAEWbW38x2B/4J/BeYFnGMF8zs/IjXN5nZgWbWycz+CDwN/AKs394jIoGKDBzP\nPAOrVvmwcdBB8NJLkMhHWaTznDQisk4sLRp742+DlOJ/78cBZcBoAOfcjcCtwF340SabA0c4536O\nOEYusHXE6+2BR4HFwGPAV8D+zrlvYiifiDSSmZ8ZtrQUpk2DNWt8/42DDoIXX0xM4MiEOWlERHOd\niEg9OAfPPuv7cJSWQs+evg/HYYf5UJKYMqjjp0g8aK4TEUk5ZnDkkfCf/8C//gU//eSfwdGrF8yZ\nk5gWDoUMkfSkoCEi9WYGf/4zvPWWb+H4+Wf/lNEDDoB//zuxfThEJD0oaIhIg5mtmxn2uef8NPV9\n+/rHm8+evenAkY63bEUkNgoaIhIzs3Uzw86c6QNGv35+ArdZs9YPHJogTaRpUtAQkUYz8wGjpMQH\njJoAsv/+PoCsXq0J0kSaKgUNEYkbM38L5Y03/C2U5s39LZa8nb+h56Ic+mqCNJEmR0FDROLOzHcS\nff11Pypl7fdfMdL9i315i2f503oTH2mCNJHMpqAhIoExg8MOcxz4+6N5nsNoyU8cybPsy1vM4Mjf\n5jKpmSBNRDKPgoaIBMrM+GGzLA7lRV7lQF7gUDbnRwYwg715m2n0Z01zTZAmkqkUNEQkcDUTpBlw\nKC/xCgfxIoeQTSUDmc6iH15l2jQ9h0MkEyloiEjgak+QZsDBvMylocPokXMaO+dty8CB0KOHn8xN\ngUMkcyhoiEjgNjZB2kvvTuSVV5rz8suw5ZZw9NGw115+uvrq6mSXXEQaS5OqiUjCbWyCtNdeg9Gj\n4YUXYI89/ORtAwf6KexFpPE0qZqIZLyNdfzs1Quef94HjnbtYNAg6N4dnnxSLRwi6UhBQyRBYmk9\nTPUWxyDL17OnfwbH669D+/Zw3HGw557wxBMKHCLpREFDJECxzO+R6nOCJLp8NTPDzp0L220Hf/mL\nv6Xy+OMKHCJpwTmXdl9APuBKS0udSKpavXq1O7xbNzczFHLVfiCFqwY3MxRyh3fr5lavXh2XfRIp\nFcr3xhvO9e3rHDjXtatzjz3m3C+/BH5akYxRWlrqAAfkuwR8ZqtFQyQgY4cPZ2h5Of0aML9HLPsk\nUiqUr2Zm2JIS2GEHOOEE2H13eOwxP129iKQWBQ2RgMydMYO+dbTt1zW/Ryz7JFIqla9mZtg334Sc\nHCgshN12g0cfVeAQSSUKGiIBcM7RuqqKusZWRJvfI5Z9EilVy7fffvDcczBvHnTuDCeeCN26wSOP\nKHCIpAIFDZEAmBlrs7Ko6yPXAWuz1p/fI5Z9EinVy7fvvvDss/DWW7DTTnDSSdC1Kzz8MPzyS1KK\nlFKSFVBFFDREAlIzv0c0s0Iheg4YEJd9EinVywewzz4wYwb85z+w887wt7/5wPHQQ00vcKT6CCZp\nIhLR4zTeX2jUiaSBmhEaz9UaofFcPUadNGSfREr18kXz9tvODRjgR6l06eLc5MnOVVUlu1TBS4UR\nQpKaNOpEJENsbH6PqSUlZGdnx2WfREr18kXTowdMmwalpb7vximnQF4ePPhgZrdwpMIIIRHQXCci\nCeNc3fN7xHOfREr18kUzfz5cdZWftK1zZxgxwvfnyMpKdsniq3duLnMqKqJ23nVAn5wc5ixfnuhi\nSQrQXCciGSqWD+RU/xBP9fJF0707PPWUDxzdu8Ppp8Muu8D990NVVbJLFx8uRUcISdOkoCEiTdKe\ne8LUqbBgAeTnwxln+MBx333pHzhSfYSQNC0KGiKSserzF/see/iZYWsCx5ln+tEq99wDP/+cgEIG\nJB1GCEnToKAhIhkl1iGdNYHj3Xf9ENmzz/aB4+670zNwXHzNNdycl8fMUOi3lg0HzAyFGJ+Xx7Ax\nY5JZPGlCFDREJGNUVlYyqKCAgkmTmFNRwbQVK5hTUUHBpEkMKiio1/Mjdt/dzwz73nv+qaPnnusf\nAHbXXekVONJxhJBkJo06EZGMMbK4mIJJk+gXZT6WmaEQ84qKGDVhQoOO+f77MGYM/N//wfbbwxVX\nwGmnQYsW8Sp1YqTjCCEJhkadiIjEKIhJ37p1gylTYOFC6NkTzj8funSBO+6A//2vsSVOHIUMSRYF\nDRHJCEEP6eza1c8M+/77cOCBcMEFPnDcfnt6BQ6RRFPQEJGMkKghnXl5fmbYRYvgoIPgwgthxx1h\n0iT46adGHVokIyloiEjGSOSQzl139TPDLloEhxwCxcU+cNx2W/0DRzr2kRNpKAUNEckYyRjSucsu\nfmbY8nLo3RsGD/aBY+JE+PHHDbfXjKrS1ChoiEjGSOaQzp13hsmTYfFiOPxwGDrUB44JE9YFjngM\nvxVJNxreKiIZK5lDOpcsgWuu8a0df/gDXHopfPnhMA6865a4Dr8VaSgNbxURiZNkDuns0gUeeAA+\n+ACOOAKCwGK3AAAYTUlEQVQuvhjG3XUZ71cP5gc232D7WIffiqQ6BQ0RkQDtuKOfGfaDDxztW77A\npdxILssZx1DW0uq37TSjqmQqBQ0RkQTYcUdjp3aX8yE7MYDpXMb15LKcsQxjLa00o6pkLAUNEZEE\nOaB/fz4MfcI9nM1H7MRRTONyriOX5Zxhf2fffscmu4gicaegISKSIJHDbzvxMfdwNh+yE3vzDA8y\nhrufuIEbboA1a5JdUpH4UdAQEUmQaMNvz84x9i1+n/fe+x/HHhviyishJweuvx402lUygYa3iogk\nSbTht5984kPGvffCFlvAsGFQVAQ1jwDRLKzSWBreKiLSREQLDDvs4CdqW7oUjj8eRo2CTp0chxbM\n4OBOu+tpopJ2FDRERFJQx45+orYFC9bQxqbw2pt9ePeTV9h3xSlMrfhGTxOVtKGgISKSwqbcfgV3\nf/83lrMjJ/EwV3MluVTwdvXlnL1oBeNGjEh2EUU2SkFDRCSFzZ0xg77V1WzPCiYymKXhwDGGEZzl\nlvHQPzuxalWySylSNwUNEZEU5ZyjdVUVkT05OvAZExnMMjpzCpOp+P58cnIco0fD998nragidVLQ\nEBFJUWbG2qwsoo0N3I7PGc8Q/rj9QZxyinH99X5Y7KhRChySWhocNMysl5lNN7MVZlZtZgOibHOV\nmX1mZj+Y2Rwz61KP415gZsvN7Ecze9PM9mlo2UREGivVhvwf0L8/s0PRL9WzQiEOO2Z/brkFli2D\n00+HG26ATp3gH/+A775LcGFFooilRaM1MB+4ADYM2mZ2KVAEnAPsC6wFZpvZZnUd0MyOB8YBI4G9\ngAXhfbaOoXwiIg1SWVnJyOJieufmptzw0cinidZccB1+WvnxeXkMGzMGgG23hZtvhuXL4cwzYexY\n38Jx5ZXw7bfJKr1IIx/YZWbVwEDn3PSIZZ8BNznnxodfbwGsBE5xzj1ex3HeBOY55waHXxvwKTDR\nOXdjlO31wC4RiYvKykoGFRQwtLycvtXVGP6DfHYoxM15eUwtKSG75mlZSSzjuBEjmDt9Oq2qqvgh\nK4sDBgxg2JgxdZZt5Uq46Sb/TI7mzaG4GIYOhd/9rv7n1cPBMlNaP7DLzHKB9sALNcucc6uBeUBB\nHftkAT1q7eOA5+vaR0QkXsYOH87Q8nL6hUMG+Cnb+1VXM6S8PCWGj2ZnZzNqwgTmLF/OM59+ypzl\nyxk1YcJGA9A22/hWjeXL4ZxzfGtHTg4MHw7ffFP3uVK5dUfSU7w7g7bH/zGwstbyleF10WwNNGvg\nPiIicVEzfDSaftXVzJ0+Peq6ZGloC8M22/iWjYoKOPdcuOUWHziuuAK+/nr9bWtadwomTWJORQXT\nVqxgTkWFHg4mjdI8QeepaY2M6z5Dhgyhbdu26y0rLCyksLCwgacSkaYo2vDRSAa0qqrKiFsI7drB\njTfCJZfAuHEwcSLcequfR2XYMNh66/Vbd2rUtO64cOvOqAkT4l62TKjfVDVlyhSmTJmy3rJVCX7w\nSlz7aIRvnSwFujvn3o3Y7mXgHefckCjHyAJ+AAbV6uvxINDWOXd0lH3UR0NE4qJ3bi5zKiqihg0H\nHJ6Tw/PLlye6WIH7+msfOG691b8uKoI3Hu3BK5+W1VkXfXJymBOnuqisrGTs8OHMnTGD1lVVrM3K\n4oD+/bn4mmuS3icm06V1Hw3n3HLgC+CwmmXhzqD7AW/UsU8VUFprHwu/jrqPiEi8bGr4aM8BG4zg\njyrVhsVuytZbw3XX+VsqxcUwaZLjjf++ymVcz1dsOOAvsnWnsXSLpmmJ5Tkarc1sTzPrHl7UOfy6\nY/j1LcAIM+tvZrsD/wT+C0yLOMYLZnZ+xGFvBs42s5PNbFfgTqAV8GDD35KISP3Vd/hoNJnQcXLr\nreHaa6Giwui4xYPczvnkUMEl3MiX/OG37RywNisrLrc40qEDrsRPLC0aewPv4FshHP75F2XAaIDw\ncNRbgbvwo002B45wzv0ccYxcWBeZw8NehwFXhY+9B9DXOfdVDOUTEam37OxsppaUMK+oiD45ORzV\noQN9cnKYV1S00aGtmfZX+e9/Dyef/AH3WWeGMJ67OIdclnMxN7GSdg1q3dmUdOuAK43TqD4ayaI+\nGiISlPp2TBxZXEzBpEnrdZysMTMUYl5RUSAdJ4NUE56GlJezb3VbJnARExjM/8him98/wQtvDqJL\nlzYbPcam6s85x8COHZm2YkWd2xzVoQPPfPqpOogGJK37aIiIpLv6frhl4l/lka07J+S0ZUGHu9mr\n44Hst++rfF91Mrvv3oYhQ+Dzz9ffryG3kDY2fwvE9xaNpAYFDRGRBmrIsNh0U/vhYC9/soBX5vWj\nosK49FJ44AHo3BkuusgHjlhuIcWrA66kBwUNEZEGaip/lUeWf6ut/MywFRVw2WXw4IM+cPTuVcop\ni75vUMfOxnTAlfSjoCEiEoOm+lf5llvCyJE+cFx+OZS9txdnuCUUM4EVbLfetnXdQoq1A66kJ3UG\nFRGJQWTHyX4Rk7HNCv9V3hQ+MJ1z/LlDVwo+P5abGcqPbM5Z3MNlXE8HPgPq17FTTwZNLHUGFRFJ\nA/qr3N9a+bnFT4xgDBXkMIIxPMKJdGYZRdzKp3So1y0khYzMphYNEZE4aKp/ldce5ruabG6jiHEM\nYzVt6L77Wzz1bC86dtzEgRKkqf4/RVKLhohIGmqqH161O3ZuQSWXcx13W2dy2t3BshUH0KULnHce\nfPJJcsqYCU9wTWcKGiIiErO6biG9d+GplC05g4qKEKNHwxNPQJcufqr6jz9OXPni8QTXdGz5TyW6\ndSIiInFT162JNWtg0iQYOxZWrYLTToMrroBOneJ/rkixPsE1k2eX1a0TERFJW3V98LdpA5deCsuX\nw5gx8NRTvoXjrLP8svpq6G2QWJ7gmmnz2CSbgoaIiCRMmzbw97/753Bcdx1Mnw477wxnnrnpwNHQ\nABDrE1w1u2x8KWiIiEjCtW4NF18My5bB9dfDjBk+cJxxhl8WTUMDQKxPcM3EeWySSUFDRESSpnVr\nGDbMt2bccAM8+6wPHKefDkuXrr9tLAGgoU9wzeR5bJJFQUNERJKuVSsYOtS3Ztx0Ezz3HOyyi+80\numRJ7AGgofOqNJV5bBJJQUNERFJGq1YwZIgPHGPHwqxZsOuucNppxld0aXAAiOUJrk11HpugaHir\niIikrB9/hLvv9rdVvviimkPdw9zB1ezEkvW229hQ1Uj1GRKb6fPYaHiriIhI2Oabw+DBvr/GDTf8\nzGvN+7ErizmZyXzITg2eXr4+tzw0j018qUVDRETSxldfVXLKcbN58fUD+N+v7WjfegaDjnmf6yYV\nBxYAMm1+FLVoiIiI1OEPf8jmuZeP5fs123LrrSGabTmQOx4ZzrnnZrN4cTDnzKSQkQwKGiIiknZa\ntoSiImPpUrj1Vnj1VejaFf76VygvT3bpJJKChoiIpK0WLeD88/0Q2Ntvh9dfh27doLAQFi1KdukE\nFDRERCQDtGjhZ4b96CO44w544w3YbTc44QQFjmRT0BARkYzRogWcc44PHHfeCW++6QPH8cfDwoXJ\nLl3TpKAhIiIZZ7PN4Oyz4cMP/XM45s2D3XeHv/wF3nsv2aVrWhQ0REQkY222mZ8Z9sMP4Z574D//\ngT32gGOPhXffTXbpmgYFDRERyXiRgeO++6CsDPbcEwYNggULkl26zKagISIiTUZWlp8Z9oMP4P77\nYf586N4djjnGfy/xp6AhIiJNTlaWnxl28WJ44AF/G2WvveDoo+Gdd+JzjnR88nYQFDRERKTJysqC\nU0/1gePBB/3IlPx8OOoof3uloSorKxlZXEzv3FwGduxI79xcRhYXU1lZGe+ipw0FDRERafKaN4dT\nTvFPFZ082f/bowcMGAClpfU7Rs2srwWTJjGnooJpK1Ywp6KCgkmTGFRQ0GTDhoKGiIhIWPPmcPLJ\n/iFfDz3k+3LsvTf07w9vv73xfccOH87QiKnlAQzoV13NkPJyxo0YEXTxU5KChoiISC3Nm8NJJ/nA\n8fDDfrTKPvvAkUf6IbLRzJ0xg77V1VHX9auuZu706QGWOHUpaIiIiNShWTM48UQfOB55xM+psu++\n8Oc/w1tvrdvOOUfrqirqmufVgFZVVU2yg6iChoiIyCY0a+Znhn3/fXj0UVi2DPbbD/70J//UUTNj\nbVYWdcUIB6zNymqSU84raIiIiNRTs2Z+ZtiFC2HKFPj4Y9h/f+jXD3bY+0Jmh6J/rM4Kheg5YECC\nS5saFDREREQaqFkzPzPse+/B//0ffPopPPDkUE7c/DXG2QG/tWw4YGYoxPi8PIaNGZPMIieNgoaI\niEiMQqF1E7U9/jhss8N+XOxeZ+uWr9Fz66Pok5PDvKIippaUkJ2dneziJoWChoiISCOFQnDccbBw\nYTOeeAK269KTuV8/g9txOb2Pm9BkQwYoaIiIiMRNKORnhl2wAJ58Er78Enr1gsMOg9deS3bpkkNB\nQ0REJM5CIT8z7Pz5MHUqfP01HHggHHoovPJKskuXWAoaIiIiAQmF/Myw77wDTz0F334LBx8MhxwC\nL7+c7NIlhoKGiIhIwEIhPzNsWRk8/TR8/70PGwcfDC+9BJn8HC8FDRERkQQJhWDgQB84nnkGVq/2\nt1MOPhhefDEzA4eChoiISIKZ+anoS0th+nRYs8Z3GL3kkmSXLP6aJ7sAIiIiTZWZnxn2yCPh2Weh\nfftklyj+FDRERESSzMyHjUykWyciIiISGAUNERERCYyChoiIiARGQSMDTJkyJdlFSAmqh3VUF57q\nYR3Vhad6SLxAgoaZtTGzW8yswsx+MLPXzWzvjWx/kJlV1/r61czaBVG+TKNfHE/1sI7qwlM9rKO6\n8FQPiRfUqJP7gK7AicDnwN+A580szzn3eR37OGBnoPK3Bc59GVD5REREJAHi3qJhZi2BY4BLnHNz\nnXPLnHOjgSXAeZvY/Svn3Jc1X/Eum4iIiCRWELdOmgPNgP/VWv4j0HMj+xkw38w+M7N/m9kfAyib\niIiIJFDcb50459aYWQlwpZktBlYCfwUKgI/q2O1z4BzgbaAFcBbwspnt65ybH2X7lgDl5eXxLn5a\nWrVqFWVlZckuRtKpHtZRXXiqh3VUF57qYb3PzpaJOJ+5AGZwMbNc4H7gIOAXoAz4EMh3zu1Wz2O8\nDHzsnDslyrq/Ao/ErcAiIiJNz4nOuUeDPkkgnUGdc8uBQ8xsc2AL59xKM3sMWN6Aw7wFHFDHutn4\njqYVwE+NKauIiEgT0xLIwX+WBi6QFo0NTmK2FbAMuNg5d1899/k3sNo5d2yghRMREZHABNKiYWZ9\n8J07PwB2Am4EyoEHw+uvBTrU3BYxs8H41o738UnrLOAQ4PAgyiciIiKJEdRzNNoC1wEdgG+BJ4ER\nzrlfw+u3BTpGbL8ZMA7YDvgBeBc4zDn3akDlExERkQRIyK0TERERaZo014mIiIgERkFDREREApO0\noGFmvcxsupmtCE+iNqDW+tZmdpuZfRqemO19Mzun1jbbmNlDZva5ma0xs1IzO6bWNluZ2SNmtsrM\nvjOze82sdSLeY33Voy7amdmD4fVrzew5M+tSa5sWZjbJzL42s0oze7L2pHRm1tHMng0f4wszu9HM\nUiZsNrYewv/XE81scXj9x2Y2wcy2qHWclK4HiM/PRK3tZ9ZxnJSui3jVg5kVmNkL4evEKjN72cxa\nRKxP6etEnK4RmXK9vNzM3jKz1Wa20syeNrOda20Tl+uhmR0crqefzOxDM9vguU7JEo96MLM9zOxR\nM/vE1n3OFkc5V6PqIZkXlNbAfOAC/IRqtY0H+uCfKrorcAtwm5kdGbHNQ/hRLUcCuwFPAY+b2Z4R\n2zwK5AGHAX8GDgTuius7abxN1cU0/Jjn/kB34BP8JHWbR2xzC/79DcK/x+2AqTUrw79Az+E7AO8P\nnAKcClwV13fSOI2th+3wHY2H4n8eTgH6AffWHCBN6gHi8zMBgJkNAX6tfZw0qYtG14OZFQAzgVnA\n3uGv24DqiOOk+nUiHj8PmXK97AXcCuwH9AaygH/H+3poZjnAv4AXgD2BCcC9ZpYqoyFjrYenItb3\nAL7EP5eqK3ANcJ2ZnV+zQVzqwTmX9C/8L/yAWsveA4bXWvY2cFXE60r8k80it/kaOD38fV742HtF\nrO+Lf1pp+2S/7/rUBf7CUA3sGrHM8I92r3mfW+Dnljk6YptdwvvtG359BFAFbB2xzTnAd0DzZL/v\neNRDHcc5Fj/PTigd66GxdRG+MHwMtItynLSqi1jrASgBRm3kuLum03WiEfWQcdfLcBm3Dpe7Z/h1\nXK6HwA3Au7XONQV4LtnvOV71UMdxbgOej3jd6HpImSbSKN4ABpjZdgBmdgj+FyrySWZzgePDzX1m\nZifg50p5Obx+f+A759w7Efs8j/+LYL+Ayx8vLfDl/W2SOuf/p//Huknq9sYn8xcitvkA/1dNQXjR\n/sB7zrmvI449Gz8UuVtQhY+j+tRDNFviH/xW89drutcD1LMuwn/ZPApc4KLPhpzudbHJejCzP+B/\n1782s7nhJvKXzSzyqcMFpPd1or6/G5l6vdwSX8Zvw697EJ/r4f7490+tbQpITbHUQzRtI44BcaiH\nVA4aF+If8vVfM/sZ38x1gXNubsQ2x+OfwfEN/pfqDnx6WxZe3x7fLPQb55/l8W14XTpYjP/BuM7M\ntjSzzczsUmB7/G0CgG2An51zq2vtu5J177N9+HXt9ZAedVGfeliPmW0NjGD9pt90rweof12MB153\nzv2rjuOke13Upx46h/8dif856Iufe+kFM9sxvC7drxP1/XnIuOulmRn+9sDrzrlF4cXtic/1sK5t\ntrCI/j2poBH1UPs4fwT+Qv2umfWuh1QOGsX4FH0kkA8MA243s0MjthmDT1+H4tPbzcATZrapv8aM\n6Pc5U45z7hfgGGBn/C/8Gvxkdc/h77tvTH3fZ8rXRUPrwcyygWeBhcDo+p4mLoUNWH3qwnxnwUOB\nIbGepvElDVY9fyZqrnF3Ouf+6Zxb4Jwbin9q8embOEVaXCca8LuRidfL2/F9CwrrsW08rodWj22S\nodH1YGa7Ac/gbzO+sMFeGx6DaMeJJqgngzaKmbXEd0o5yjk3K7x4oZntBVwMvGhmnfEdo7o65xaH\nt3nPzA4MLz8f+AJ/bzry2M2ArdgwoaWscFNmfvjDczPn3Ddm9ibwn/AmXwCbmdkWtdJrO9a9zy+A\nfWodepvwv2lRF/WoBwDMrA2+ae974Bi37om0kAH1APWqi0Pwf82v8n/s/OYpM3vVOXcoGVAX9aiH\nz8P/ltfatRzYIfx92l8nNlUPmXi9NLPbgD8BvZxzn0Wsauz18IuIf7eptU07/K3Ynxtb/nhpZD3U\nHKMr/vbInc6562qdotH1kKotGlnhr9pp6VfWlblVeP3GtikBtgwHlBqH4dPYvHgWOBGcc5XhC8hO\n+H4Zz4RXleI7bB1Ws214mNMO+L4u4Oti9/DthBp9gFXAItLIRuqhpiXj3/gOoAOi/CJkTD3ARuvi\nOmAPfGfQmi+AwcBp4e8zpi7qqgfnXAXwGb4TXKSd8Z1kIYOuExv5ecio62X4w/Uo4BDn3Ce1Vjf2\nelgesc1hrK9PeHlKaEQ9lEQs6wa8CDzgnPtHlNM0vh6C6gG7qS/8cK098UOxqoGLwq87hte/hJ/z\n5CD8sK1T8fOgnB1e3xz4EN+RaR/8X2/DwhXbN+I8z+FHq+yDn3b+A+ChZL3vGOvi2HA95OJ/qJYD\nj9c6xu3h5Qfjm0XnAq9FrA8BC/DD/PbA36teCVyd7Pcfr3oA2gBv4ocB5uJTeM1XzaiTlK+HeP1M\nRDlm7dEKKV8XcfrdGIwfTTAI2BG4GlgL5EZsk9LXiTj8bmTS9fL28P9nr1q/4y1rbdOo6yH+c2cN\nftTFLvhWn5+B3smugzjWQzd8v5x/1jpG5GicRtdDMivpoPAvzK+1vu4Pr28H3Ad8Gr4oLAIG1zrG\njsAT+ObRSuAd4K+1ttkSeBifVL8D7gFaJfuHpIF1cSG+s9dP4R+aUdQafojvPX4rfrhaZbhe2tXa\npiN+PPSa8C/VDYQ/gFPhq7H1EN6/9r41x9shXeohXj8TUY75KxsOI0/puohXPQB/x7dgVAKvAwW1\n1qf0dSJO14hMuV5Gq4dfgZMjtonL9TBc76X4FtKPgL8l+/3Hsx7wnaSjHWNZPOtBk6qJiIhIYFK1\nj4aIiIhkAAUNERERCYyChoiIiARGQUNEREQCo6AhIiIigVHQEBERkcAoaIiIiEhgFDREREQkMAoa\nIiIiEhgFDREREQmMgoaIiIgE5v8BW0/uc689+/YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x104989550>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(x,t,'ro')\n",
    "plt.plot(testx,testt,'b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Phew again -- it looks just like the old version. So, what's the point?\n",
    "\n",
    "Well, to add, say, an $x^2$ term we just need to add a column to $\\mathbf{X}$ (and the test version):\n",
    "\n",
    "$$ \\mathbf{x}_n = \\left[\\begin{array}{c} 1 \\\\ x_n \\\\ x_n^2 \\end{array}\\right] $$\n",
    "\n",
    "$$ \\mathbf{X} = \\left[ \\begin{array}{c} \\mathbf{x}_1 \\\\ \\mathbf{x}_2 \\\\ \\vdots \\\\ \\mathbf{x}_N\\end{array}\\right] = \\left[ \\begin{array}{ccc} 1 & x_1 & x_1^2 \\\\ 1 & x_2 & x_2^2 \\\\ \\vdots & \\vdots & \\vdots \\\\ 1 & x_N & x_N^2 \\end{array}\\right] $$\n",
    "\n",
    "The equation for $\\mathbf{w}$ doesn't change **at all**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  4.55597856e+02]\n",
      " [ -4.43160485e-01]\n",
      " [  1.10151552e-04]]\n"
     ]
    }
   ],
   "source": [
    "X = np.hstack((np.ones_like(x),x,x**2))\n",
    "XX = np.dot(X.T,X)\n",
    "invXX = np.linalg.inv(XX)\n",
    "Xt = np.dot(X.T,t)\n",
    "w = np.dot(invXX,Xt)\n",
    "print w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "testx = np.linspace(1896,2012,100)[:,None]\n",
    "testX = np.hstack((np.ones_like(testx),testx,testx**2))\n",
    "testt = np.dot(testX,w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1063e4c50>]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhoAAAFkCAYAAABmeZIKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3X2c1XP+//HH69QUZbbCCqlmIkwRJtIgoUtWhbBi136t\nXX7Ujo2wu/VVbK1lVbKy6+vaxlgrq1rSpiQyhRnXJnYpF0lyVVNt2+i8f3+8z+h0mqmZM+dzLp/3\n2+1za87n8nXezXzO67w/7wtzziEiIiIShFCqAxAREZHspURDREREAqNEQ0RERAKjRENEREQCo0RD\nREREAqNEQ0RERAKjRENEREQCo0RDREREAqNEQ0RERAKjRENEREQC0+hEw8z6mNlsM1tlZmEzGxq1\nrbmZ3WRmb5jZhsg+D5jZfg0470gzW2Fm/zGzpWZ2TGNjExERkfQST41Ga+A1YCQQO1FKK+BI4Hrg\nKOBM4BBg1s5OaGY/BCYD4yPHvQ7MM7O944hPRERE0oQ1ZVI1MwsDZzjnZu9kn6OBZUBn59wn9eyz\nFFjmnLsi8tqAj4HbnHM3xx2giIiIpFQy2mi0xdd8fFPXRjPLA3oCC2rXOZ/9PAOUJCE+ERERCUjz\nIE9uZi2B3wMPO+c21LPb3kAzYE3M+jX4xy51nXcvYBCwEtickGBFRERyw25AATDPOfdl0BcLLNEw\ns+bA3/C1GZfHcwp2bANSaxDwUJyhiYiICFwAPBz0RQJJNKKSjI7AKTupzQD4AtgKtI9Zvw871nLU\nWgkwY8YMioqKmhZsiv2/00/nT6tXY3Vsc8Bl++3Hn//xj52eY/To0UydOjWQ+DKJymEblYWncthG\nZeGpHKCqqoof/ehHEPksDVrCE42oJKMLcLJz7uud7e+cqzGzCqAfMDtyDou8vq2ewzYDFBUVUVxc\nnKjQU+LUs89m7fTpDA6Hd9g2NxTitHPO2eV7bNOmTcaXQyKoHLZRWXgqh21UFp7KYTtJaXoQzzga\nrc3sCDM7MrKqS+R1RzNrBswEioEfAXlm1j6y5EWdY4GZRT9OmQJcYmYXmtmhwJ/xXWXvj/N9ZYwx\nkyYxpaiIuaHQd8+JHD7JmFpUxFUTJ6YyPBERkSaJp0bjaOBZ/Oehw49/AfAAfvyMIZH1r0XW17a1\nOBlYHFlXiG8ECoBz7tHImBk34B+hvAYMcs6tjSO+jJKfn8/M8nImjxvHlNmzaVVTw6a8PI4fOpSZ\nEyeSn5+f6hBFRETi1uhEwzn3HDuvCdllLYlzrksd6+4A7mhsPNkgPz+fCdOmwbRpOOfwT45EREQy\nn+Y6STPxJBkjRowIIJLMo3LYRmXhqRy2UVl4Kofka9LIoKliZsVARUVFhRr1iIiINEJlZSU9e/YE\n6Omcqwz6eqrREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo\n0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjR\nEBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQ\nERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAR\nEZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBER\nkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGR\nwCjREBERkcA0OtEwsz5mNtvMVplZ2MyGxmw/08yeNrO1ke09GnDOn0T23Rr5N2xmmxobm4iIiKSX\neGo0WgOvASMBV8/2F4Br69len3XAvlFL5zhiExERkTTSvLEHOOeeBp4GMDOrY/uMyLbOwA7bd35q\nt7ax8YiIiEj6Sqc2GnuY2Uoz+8jMnjCzbqkOSERERJomXRKNd4GfAkOBC/BxvWhmHVIalYiIiDRJ\nox+dBME5txRYWvvazMqBKuASYHx9x40ePZo2bdpst27EiBGMGDEioEhFREQyR1lZGWVlZdutW7du\nXVJjMOca014z5mCzMHCGc252Hds6AyuAI51zb8Rx7keBGufcBXVsKwYqKioqKC4ujiNyERGR3FRZ\nWUnPnj0BejrnKoO+XtCPTuLKYswsBBwGrE5sOCIiIpJMjX50YmatgYPY1qOki5kdAXzlnPvYzNoB\nnYAOkX0OjfRO+cw5tyZyjgeAVc6530Re/y/+0cm/gbbANfjurXc35c2JiIhIasXTRuNo4Fl8bYUD\nJkfWP8C2Bp33RW2vfTh0PXBD5OeOwNaoc7YD/g8/fsbXQAVQ4pxbHkd8IiIikibiGUfjOXbyyMU5\n9wA+6djZOU6JeX0lcGVjYxEREZH0li7dW0VERCQLKdEQERGRwCjREBERkcAo0RAREZHAKNEQERGR\nwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREZHA\nKNFIE1u3wk03wYYNqY5EREQkcZRopIn33oOJE2HQIPjmm1RHIyIikhhKNNJEUREsWABVVXDyybB2\nbaojEhERaTolGmmkVy947jlYvRpOPBE++STVEYmIiDSNEo00c/jh8PzzsGkT9OkD77+f6ohERETi\np0QjDXXt6pONFi3ghBPgrbdSHZGIiEh8lGikqU6dYPFiaN/eP0ZZtizVEYmIiDSeEo001r49LFrk\nG4r26wcLF6Y6IhERkcZRopHm2raFf/4Tjj8eTjsNnngi1RGJiIg0nBKNDNC6NcyeDUOHwvDhcP/9\nqY5IRESkYZRoZIiWLaGsDC6+GC66CKZMSXVEIiIiu9Y81QFIwzVrBnfeCXvtBVddBV9+6UcTNUt1\nZCIiInVTopFhzODGG2HPPeGaa3yyMX26T0JERETSjRKNDHX11bD33vCzn/lkY8YM/3hFREQknaiN\nRga76CJ4/HGYM8f3SKmuTsx5nXOJOZGIiOQ8JRoZbtgw3/31lVfgpJPg88/jO091dTXjS0vpX1jI\nGR070r+wkPGlpVQnKnsREZGcpEQjC5x4op+MbdUqP97GihWNO766uprhJSWUTJ/O/JUrmbVqFfNX\nrqRk+nSGl5Qo2RARkbgp0cgSRx4JL77ofz7uOHj99YYfe8vYsVxZVcXgcJjaDiwGDA6HGV1VxeRx\n4xIdroiI5AglGlmkSxdYsgQ6dPC1HIsWNey4JXPmMCgcrnPb4HCYJbNnJy5IERHJKUo0ssw++8Cz\nz0KvXjBoEDz22M73d87RuqaG+obiMKBVTY0aiIqISFyUaGSh/Hx48kk/XPm558Ltt9e/r5mxMS+P\n+tIIB2zMy8M0KpiIiMRBiUaWatHCj60xejT84hfwm99AfZUSxw8ZwrxQ3b8KT4dCnDB0aICRiohI\nNtOAXVksFILJk2H//WHMGPj0U7jrLsjL236/MZMmMXzhQlxUg1CHTzKmFhUxc+LEVIQvIiJZQDUa\nOeCqq+Chh+Dhh+H003cc2Cs/P5+Z5eUsGzWKgQUFDOvQgYEFBSwbNYqZ5eXk5+enJnAREcl4qtHI\nEeefD/vuC2eeCX37wlNP+de18vPzmTBtGkybhnNObTJERCQhVKORQ045BZ5/HtasgZISePfduvdT\nkiEiIomiRCPH9OgB5eXQqpUf2GvJklRHJCIi2UyJRg7q1AleeAEOPxz69dv1WBsiIiLxUqKRo9q1\ng3nzfJuNc8+FqVNTHZGIiGQjNQbNYS1b+t4onTvDlVfCypUwZQo0a5bqyEREJFso0chxoRD8/vc+\n2Rg1Cj780HeDbdUq1ZGJiEg20KMTAeCyy2D2bHjmGTjpJN8zRUREpKmUaMh3fvADWLwYPv4YeveG\nqqpURyQiIplOiYZsp7gYli2D1q39WBsLF6Y6IhERyWRKNGQHnTr58TVqp5q/775URyQiIplKiYbU\nqU0bP9X8RRfBT38KY8dCOJzqqEREJNOo14nUKy8P7rwTunaFa66Bf/8b7r8fdt891ZGJiEimUI2G\n7JQZXH21Hz10zhw4+WT1SBERkYZrdKJhZn3MbLaZrTKzsJkNjdl+ppk9bWZrI9t7NPC855hZlZn9\nx8xeN7NTGxtbrnLOBX7M8OHw3HN+nI1jj4W33mr0JUVEJAfFU6PRGngNGAnU9WnVGngBuLae7Tsw\nsxLgYeAu4EjgCeAJM+sWR3w5obq6mvGlpfQvLOSMjh3pX1jI+NJSqqurE3pMtGOOgZde8u03jjsO\n5s5N1LsREZGs5ZyLewHCwNB6tnWObO/RgPM8AsyOWVcO3FHP/sWAq6iocLlo/fr1bkD37m5uKOTC\n4By4MLi5oZAb0L27W79+fUKOqf/6zp1+unOhkHO33upcOJzIdyciIkGqqKhw+IqAYteEHKChS7q0\n0SgBnolZNy+yXmLcMnYsV1ZVMTgcxiLrDBgcDjO6qorJ48Yl5Jj65OfDE0/A6NHwy1/C5ZdDTU1T\n35WIiGSjdEk09gVimxiuiayXGEvmzGFQPX1NB4fDLJk9OyHH7EyzZnDLLXD33X459VT4+utGnUJE\nRHJAOndvNXbRxmP06NG0adNmu3UjRoxgxIgRQcaVUs45WtfUfFcrEcuAVjU1OOcws7iPaaiLL4YD\nD/SNRY891vdMOeSQRp1CREQCUlZWRllZ2Xbr1q1bl9QY0iXR+AxoH7NuH3as5djO1KlTKS4uDiyo\ndGRmbMzLw0GdiYMDNublbZcwxHNMY5x0km8kOmSITzYefRQGDozrVCIikkB1ffmurKykZ8+eSYsh\n6EcnDe1DWQ70i1k3ILJeYhw/ZAjzQnX/1z0dCnHC0KE7rI/nmMY48EAoL/e9UU49FW67zbc4FRGR\n3BbPOBqtzewIMzsysqpL5HXHyPZ2ZnYE0B3/BfrQyPb2Ued4wMx+F3XaacCpZnalmR1iZhOAnsDt\ncb6vrDZm0iSmFBUxNxT6LpNzwNxQiKlFRVw1cWJCjmmsNm38o5Nf/hKuuAIuuQS2bGnyaUVEJIPF\nU6NxNPAqUIH/rJoMVALXR7YPjWyfE9leFtl+adQ5OhLV0NM5Vw6MAC7Bj9FxFjDMOfdOHPFlvfz8\nfGaWl7Ns1CgGFhQwrEMHBhYUsGzUKGaWl5Ofn5+QY+LRrBlMnuwnYnvwQejXDz7/PCGnFhGRDGQu\nA+u3zawYqKioqMi5Nhp1iacRZzzHNFZ5OZx5JrRsCbNmwZFH7voYEREJVlQbjZ7Oucqgr5cu3Vul\nCeJJGIJOMgBKSuCVV+D734fjj/eNREVEJLco0ZBAHXAAPP88nHEG/PCH8JvfwNatqY5KRESSJV26\nt0oW2313mDEDjjoKrr0W3ngDHnrINx4VEZHsphoNSQozGDMGnnwSXngBevWC5ctTHZWIiARNiYYk\n1eDB8PLL0Ly5TzYaOfK5iIhkGCUaknRdu8LSpdC/PwwbBtdfD/VMwyIiIhlOiYakRH4+PPYYTJzo\nE40zz4QkD78vIiJJoERDUiYUgrFj/Wiizz3nH6W8oyHaRESyihINSbkf/MCPt9GihU82Hnss1RGJ\niEiiKNGQtHDQQX4k0dNPh3PO8d1gv/021VGJiEhTKdGQtLHHHlBW5udKmTzZTzWveVJERDKbEg1J\nK2Zw5ZWwYAG8/Tb07AnLlqU6KhERiZcSDUlLfftCZSV07Ah9+sAdd0AGzv8nIpLzlGhI2urQARYt\ngksvhZEj4cc/ho0bUx2ViIg0hhINaTAXZ5VCvMeB74nyxz/Cww/D3/8Oxx4L776b+Os0VjKvJSKS\nyZRoyE5VV1czvrSU/oWFnNGxI/0LCxlfWkp1dXUgx9VnxAh46SXfE+WYY+BvfwvmOjuTzGuJiGQN\n51zGLUAx4CoqKpwEZ/369W5A9+5ubijkwr6JhAuDmxsKuQHdu7v169cn9LiGxeTcuec6B85ddtl/\nXb+iIwK5TjLfk4hIMlVUVDjAAcUuCZ/ZqtGQet0ydixXVlUxOBzGIusMGBwOM7qqisnjxiX0uIbI\nz4dHHvGPU+68M8Sqqj/RLdwh4deJFeR7EhHJZko0pF5L5sxhUD2znQ0Oh1lSz9Sr8R7XUGYwahQU\ntz+bTexPMZU8xakJv060oN+TiEi2UqIhdXLO0bqm5rtv77EMaFVTs0OjyHiPiye+/UOvUEkxvVnK\nD3iKX3Ej39IsodepvVYy3pOISDZSoiF1MjM25uVR30enAzbm5WG2/cdvvMfFG9+efMVshnIT13AL\nYziZZ1nF/gm7TvS1gn5PIiLZSImG1Ov4IUOYF6r7V+TpUIgThg5N6HHxxhfCcQ1/YBEnsYJCjuQ1\nfmunJew60deqSyLfk4hI1klGi9NEL6jXSVLU9rR4KqanxVMN7HXS2OMSEd/n7OWO5ikHzv3yl5vd\nli0JuVTS3pOISNDU60TSRn5+PjPLy1k2ahQDCwoY1qEDAwsKWDZqFDPLy8nPz0/ocYmI7/yCfE77\nxTwmTtzM7be3pG9f+PDDYK4VxHsSEck25jKwAZuZFQMVFRUVFBcXpzqcnOGci6sdQrzHNfU6S5fC\neefBunVw771w5pnBXUtEJFNUVlbSs2dPgJ7Oucqgr6caDWmweD9Yk/WBHHud3r3h1Vfh5JPhrLN8\nl9jNm4O5loiI1E2JhmS1du1g5kw/++vdd/u5UpYvT3VUIiK5Q4mGZD0zuOwyP1fKli3Qsyfcc4+m\nnRcRSQYlGpIzevSAV16B88+Hn/3Mt9/45ptURyUikt2UaEhOad0a7roLHn0U5s2DI4+EJUtSHZWI\nSPZSoiE56Zxz4PXX4YAD4MQTYcIEPwW9iIgklhINyVmdO8OiRXDddfDb30LfvrBiRaqjEhHJLko0\nJKc1bw7jx8Pzz8Onn8IRR8CMGWooKiKSKEo0RIDjjoPXXoNhw+DHP/YNRr/+OtVRiYhkPiUaIhFt\n2sBf/gJlZTB3rq/dWLQo1VGJiGQ2JRoiMc47D954Aw48EE45Ba65Bv7731RHJSKSmZRoiNShUyd4\n5hm46Sa49VY45hiffIiISOMo0RCpR7NmcPXV8PLLvnHoMcfAH/4AW7emOjIRkcyhRENkF444wicb\npaVw7bV+krYPPkh1VCIimUGJhuQ814C+rLvt5msznn0WPv7YD2d+113qBisisitKNCQnVVdXM760\nlP6FhZzRsSP9CwsZX1pKdXX1To/r29e31RgxAi65BE4/3Y+/ISIidVOiITmnurqa4SUllEyfzvyV\nK5m1ahXzV66kZPp0hpeU7DLZyM/3tRn/+AdUVsJhh8HDD6t2Q0SkLko0JOfcMnYsV1ZVMTgcxiLr\nDBgcDjO6qorJ48Y16Dw/+AG89RYMGgQXXABnnw2ffx5Y2CIiGUmJhuScJXPmMCgcrnPb4HCYJbNn\nN/hce+3lB/j6299g8WLo3t3/LCIinhINySnOOVrX1HxXkxHLgFY1NQ1qIBrt7LN97UafPnDuuX5Z\nu7bJ4YqIZDwlGpJTzIyNeXnUl0Y4YGNeHmb1pSL1a98eZs6ERx6BhQuhWzfVboiIKNGQnHP8kCHM\nC9X9q/90KMQJQ4fGfW4z+OEP4e234cQTfc3G2WfDmjVxn1JEJKMp0ZCcM2bSJKYUFTE3FPquZsMB\nc0MhphYVcdXEiU2+Rvv28Nhj8Oijvu1Gt26afl5EcpMSDck5+fn5zCwvZ9moUQwsKGBYhw4MLChg\n2ahRzCwvJz8/PyHXMYNzzvG1G4MG+ennhwzxA36JiOQKa2yjt3RgZsVARUVFBcXFxakORzKccy6u\nNhmNNWsWXHYZbNgAN9/sB/yq5wmOiEhgKisr6dmzJ0BP51xl0NfTbU5yXjKSDIBhw+Cdd/w09Jdd\n5udMee+9pFxaRCRlGp1omFkfM5ttZqvMLGxmO7ScM7MbzOxTM9tkZvPN7KBdnHN85FzRyzuNjU0k\n3bVtC//3f7BgAXzyiZ8z5Xe/g5qaVEeW/jKx9lVE4qvRaA28BoyEHXsJmtm1wCjgUqAXsBGYZ2Yt\ndnHet4D2wL6R5YQ4YhPJCKecAm++6WeEve466NkTXnop1VGln3jnpBGR9NHoRMM597Rz7jrn3BNQ\n57hHVwC/dc7Ncc69BVwI7A+csYtTf+ucW+uc+zyyfNXY2EQySatWvq3Gyy9DXh707g1XXAH6DPWa\nOieNiKSHhLbRMLNCfG3Egtp1zrn1wDKgZBeHd408jnnfzGaYWcdExiaSro46CpYt89PQ33237wo7\na1aqo0q9RM1JIyKplejGoPviH6fEDk+0JrKtPkuB/wEGAf8PKAQWm1nrBMcnkpaaN4errvJdYXv0\ngDPOgLPO8u04clUi56QRkdRpnqTrGHW056jlnJsX9fItM3sJ+BA4F7ivvuNGjx5NmzZttls3YsQI\nRowY0bRoRVKkoMBPP/+3v/nHKEVFMHEijBzpk5Fc0Zg5aZLVa0gkE5WVlVFWVrbdunXr1iU1hkTf\nuj7D3wPas32txj7Aqw09iXNunZm9B+y0t8rUqVM1joZkHTM/dPnAgTB2LIweDQ88AH/+M/Tqlero\nkiN6Tpq60oimzEkjkkvq+vIdNY5GUiT00YlzbgU+2ehXu87MvgccC7zY0POY2R7AgcDqRMYnkkna\ntoXp02HpUj90ee/ecPnl8PXXqY4sOYKck0ZEkieecTRam9kRZnZkZFWXyOvaxpu3AuPMbIiZHQ48\nCHwCzIo6xwIzuzzq9R/M7EQz62xmxwF/B74Ftq/vEclBvXr5nilTp/r5Ug49FB58MPvnTUnGnDQi\nErx4ajSOxj8GqcD/3U8GKoHrAZxzNwN/BO7E9zbZHTjVObcl6hyFwN5Rrw8AHgaWA48Aa4Hezrkv\n44hPJOs0b+7bbCxfDv36wU9+AiedBG+9lerIgpOsOWlEJFia60QkAy1Y4BuI/vvfftCvCRPge99L\ndVTBUsNPkcTQXCciskv9+sEbb8CkSXDnnXDIIdk/Db2SDJHMpERDJEO1aAHXXusfp/Tp46ehP/FE\neO21VEcmIrKNEg2RDNexIzz6KMyfD19+6edNufxy/3O6ysRHtiISHyUaIlmif394/XW45RZ46CE4\n+GDfPfbbb1MdmacJ0kRykxINkSySl+cH+HrvPT+M+S9+4edSWbgwtXFpgjSR3KVEQyQLtW8P99zj\np57Pz/eNR4cPhw8+SE08miBNJHcp0RDJYkcfDUuW+Ecpy5b5uVN+9StYvz65cWiCNJHcpURDJMuZ\nwfnnw7vvwq9/DbfdBl27+inpt24N/vqNmSBNRLKPEg2RHNG6tR/Y6733YMAA+PnPffuN+fODvW70\nBGl10QRpItlNiYZIjjngAD+417Jl0KaNnyX2tNPg7beDu6YmSBPJXUo0RHJUr16weDE89piv5ejR\nAy65BFYHMGeyJkgTyV1KNERymJnvjfLOOzB5Msyc6dtvXH89bNiQuOtogjSR3KVJ1UTkO19/DTfe\nCNOmwZ57wvjxcPHFfnyORNIEaSKpo0nVRCRl2rWDm2/2PVQGDPBDmR92mK/pSOR3EiUZIrlDiYZI\nksRTe5iqGseCAnjwQXj1VejSBc4+G3r3hmefTY/4RCRzKNEQCVA883uk05wgRxwBc+fCggX+9Smn\nQL9+3/Lz8/6QFvGJSAZwzmXcAhQDrqKiwomkq/Xr17sB3bu7uaGQC/snDy4Mbm4o5AZ07+7Wr1+f\nkGOSJRx2bsaMTa5Vi/cdOHcOf3VVHJI28YlIw1RUVDh8x69il4TPbNVoiAQknvk90nlOEDN4b9m1\nPFpzKHdzMUvpTXfe5qfcy6HhjimPT0TSkxINkYDEM79Hus8JsmTOHE5zNVzMvfyLrkxlNE9xGgfz\nHnPCt7Hg8WUpjU9E0o8SDZEAuDjm94jnmGSKja8lWyjlj3xAF65nPI8wghc/WcQvf+lYsyYlIYpI\nGlKiIRKAeOb3SPc5QeqLrzWb+DW/5wMK6dz2T9x3n9GlC1xzDaxdm5JQpQ6pSlBFlGiIBCSe+T3S\nfU6QncVXHtrAhReuZMUKGD0a/vQnKCz009J/8UWSAxUgvXowSQ5LRovTRC+o14lkgNoeJE/F9CB5\nqgG9ThpzTDI1Jr4vvnDu1792bo89nGvd2rlrrnHu889TGHyOSeceTJJa6nUikiXimd8j3ecEaUx8\ne+0Fv/sdrFgBpaVwxx1+ILAxY1AbjiRI5x5Mkls014lIkjjX+Pk94jkmmRoT31dfwa23+nlUtmzx\nM8VefbWftl4Sr39hIfNXrqyzcbEDBhYUMH/FimSHJWlAc52IZKl4EoZ0TjKgcfHtuSfccAN8+CH8\n+tcwY4Yf3vySS+D99wMMMge5NO/BJLlFiYaIJFXbtnDddbByJUyaBLNmwcEHwwUXwJtvpjq67JDu\nPZgktyjREJGUyM/3j05WroTbboMXXoAePWDIEHjxxcRcI5e/sad7DybJHUo0RCSldt8dRo6Ef/8b\nHnjAP0Y5/njo0wf+8Q+oZ6DUeqlLpzdm0iSmFBUxNxT6rmbDAXNDIaYWFXHVxImpDE9yiBINEUkL\neXlw4YXw1lvwxBOwdauv3ejRwycgW7bs+hzV1dUMLymhZPp05q9cyaxVq5i/ciUl06czvKQkp5KN\ndO/BJLlDvU5EJC055x+n3HQTPPkk7L8/XHGFbzzatm3dx4wvLaVk+nQG11ENMjcUYtmoUUyYNi3g\nyNNTuvdgkuRRrxMREfxssbWPT95+GwYPhv/9X+jUCa680rftiJXuk9KlkpIMSRUlGiKS9rp1g3vu\n8cnFqFFw//1w4IFw7rlQXu73UZdOkfSkRENEMsZ++/nRRj/+GG6/HV57DY47Dnr3hkceMaqb764u\nnSJpRomGiGSc1q3hsstg+XKYPdu/Pv98eOWLl7nIfsMX7LXDMenYpVO1K5ILlGiISMYKhXzPlAUL\n4I034MzhrfgL17E/n3AR91LJUWnXpVPdbyXXKNEQkaxw+OHwwAMteP+DGk48bj5lzQbSk0ratXyZ\nuwfcS9n3pRbsAAAX90lEQVRzqe/Sqe63kouUaIhIViko2INnlgxhw+YOzJzpKD7uaB6f9xO6dctn\n3Dj46KPUxaYZVSUXKdEQkazUvDmcdZaxcCG8847voXLbbVBYCMOGwdy5flCwZFL3W8lFSjREJOsV\nFcEf/wiffgp33OFnkD3tNDjoILjxRvjss+BjUPdbyVVKNEQkZ+yxB1x6Kbz6KixdCied5Keu79gR\nzjoLnn46uFoOzagquUqJhojkHDM49li47z5fyzF1qp/U7dRToUsXmDDB13okmmZUlVykRENEclq7\ndn600ddf97UcAwfC5Mm+LcfAgfDXv8LmzYm5lmZUlVykRENEhG21HHfdBatX+yHPN22C887zE7qN\nHAmvvOIne4uXZlSVXKTZW0VEduLdd/3cKg8+6B+zdO/up7O/4ALo0KFp545nRlXNwipNpdlbRUTS\nyCGH+J4pH30ETz3lBwYbP943IB04EP7yF9iwIb5zNzRh0GiiksmUaIiINECzZr6xaFmZ7w57112+\n7caFF0L79r6GY+5c+PbbxF5Xo4lKplOiISLSSG3awMUXw+LFsGIFjB3ru8yedpqfYXbkSHjhBahn\nbK5G0WiikumUaIiINEFBAfzmN/D221BZCRdd5GeU7dPHbxszBl5+Of5GpBpNVDKdEg0RkQQwg6OO\ngptv9mNwLF4Mp5/uG5H26gUHHgi/+pVPRhqadGg0UckGSjRERBIsFPI1Gnfc4XuqzJ8P/fr5dh09\ne/qhz6+9dtc1HRpNVLJBoxMNM+tjZrPNbJWZhc1sh6HszOwGM/vUzDaZ2XwzO6gB5x1pZivM7D9m\nttTMjmlsbCIiTZXo2oHmzaF/f59kfPYZ/POf/vW99/qajs6d4Yor4Lnn6h7+XKOJSqaLp0ajNfAa\nMBJ2TLTN7FpgFHAp0AvYCMwzsxb1ndDMfghMBsYDRwGvR47ZO474REQaJVndR/PyYMAAuPNOPyjY\nggV+JtmZM/28K/vtBz/9qW/jsWmTP0ajiUqma9KAXWYWBs5wzs2OWvcp8Afn3NTI6+8Ba4CfOOce\nrec8S4FlzrkrIq8N+Bi4zTl3cx37a8AuEUmI2u6jV1ZVMSjSs8MB80IhphQVJWXEznDYP0b5+99h\n1ixYvhx2392P0zFkCPTtu4EZfxzLktmzaVVTw6a8PI4fOpSrJk4MNDYNDpadMnrALjMrBPYFFtSu\nc86tB5YBJfUckwf0jDnGAc/Ud4yISKKkQ/fRUMgPf/7730NVlU80JkyAL76ASy6Brl33YO7SafS9\neAXXzfqYf36wggnTpgWSZGhwMEm0RDcG3Rf/ZWBNzPo1kW112Rto1shjREQSIh27jx5yCFxzjR+L\nY80aeOAB35bj5pvh6KONDh38OB4zZ8I33yTuuhocTILQPEnXqa2NTOgxo0ePpk2bNtutGzFiBCNG\njGjkpUQkFzWm+2iqHiHsvbcfffTCC6GmBpYsgSef9Mu99/oRS0tKYNAgvxQX+3XxiK7dqVVbu+Mi\ntTsTpk1LzBuLokc0wSkrK6OsrGy7devWrUtqDAltoxF5dPI+cKRz7o2o/RYBrzrnRtdxjjxgEzA8\npq3H/UAb59yZdRyjNhoikhD9CwuZv3JlncmGAwYUFPDMihXJDqtBPvwQ5s2Dp5+GZ56B6mrYc0/f\nlXbgQN/wtHPnhp9vV2UxsKCA+Qkqi+rqam4ZO5Ylc+bQuqaGjXl5HD9kCGMmTdIstgHL6DYazrkV\nwGdAv9p1kcagxwIv1nNMDVARc4xFXtd5jIhIoiSq+2gqBs3q3Nm34Xj8cfjySz9I2MiRfgK4Sy/1\nI5MeeKDf569/hc8/r/9cyRwcTI9ocks842i0NrMjzOzIyKoukdcdI69vBcaZ2RAzOxx4EPgEmBV1\njgVmdnnUaacAl5jZhWZ2KPBnoBVwf+PfkohIwzWl+2g6NZzMy/ODhN1wAyxd6huSPv64nwju+efh\nvPP85G/du8OoUfDYY9snHskcHCwdGuBKEjnnGrUAfYEwsDVmuTdqnwnAp/hHIvOAg2LO8QFwXcy6\ny4GVwH+AcuDoncRQDLiKigonItJU69evd+NLS13/ggI3tEMH17+gwI0vLXXr16/f6TEDund3c0Mh\nF/YDfLowuLmhkBvQvftOj02FVauce+gh5372M+cOOsi5SMiuqMi5Sy91bsYM5674n+vc3FBo28ao\n5alQyI0vLU1ILP0KCr4rs9glDK5/QUFCriN1q6iocPjcsdg1MgeIZ2lSG41UURsNEQmKa2DDxPGl\npZRMn75dw8lac0Mhlo0aFUjDyUT55BNf07F4sR+VtKrKr2/Z/DN6f7uIs3me4yjncN7gmZBjagPH\nFNlV+TnnOKNjR2atWlXvPsM6dOCJjz9WA9GAZHQbDRGRTNfQD7d07BbbGAccACNGwJ/+BO+84x+j\nPPEEXHr5nry/77GUMpWeVLKbreey/ZZz+KAKFizIZ/XqHc/VmEdImr8l9ySre6uISNZwGdAttrG+\n/30/HPqwYS2YNq2QzZvhlVccS5e2ory8K2VlMGWK3/eAA/w8LcccA926bWLqtYO49r1lTIgeWXX6\ndIYvXFhnLcjxQ4Ywr57aIM3fkn2UaIiINFL0t/L6uoJm+rfy3XaDE04wTjhh27pPPoFly+Cll/zy\nu99BdXUr4EU+5t8UU8lRvEoxlRSHX/2uYWfsI6QxkyYxfOFCXFSDUIdPMqYWFTFT87dkFSUaIiJx\nyMVv5Qcc4Jfhw/3rcBiO79SPy1ftx6sU8ypH8Xt+xXr8QIr7hldT83/vsXl36NEDDj/cj3qan5/P\nzPJyJo8bx5SY+VtmBjx/iySfGoOKiMShdiyI0fV9K0/CZGypVlfDzjDGB3ThdY7gNY7krt2OpeU+\nA/joI1+707y5TzYOOwy6dfNLUZGja1ejRb1zfEsiJbsxqGo0RETioG/ldT9CCuE4iPc5iPc5i8cp\n39ePrPrNN/DWW3558014+21YuBDWrgUwmjXzg4sdeqhPRA4+2C9du8K++0IGP4XKeUo0RETilJ+f\n79sfTJuWUQ0/E6mhj5DatoUTTmC7Nh/gE4133vEz1i5fDu++6wcTW7nSD6wBsMce0KWLT0Rql8JC\nP/Jp586+PUlD5er/Uyop0RARSYBc/fBqasPO738f+vb1S7T//hc++ADeew/+9S94/32/PP64n+Nl\n69Zt++63H3Tq5JeOHf2/HTr49iQdOsAee1Rz63jNq5IqaqMhIiJNUl1dzeRx41gS8wjpqoAeIX37\nLaxaBStW+GXlSvj4Yz/Hy0cf+Z83b44+Ikxb1tKZ1ezHavblMzbyOa+3N66e+As6ddqdvffmu6VV\nq+2vl221IMluo6FEQ0REEiaZH8r1Xcs5+Oorn4z8ftyf2ecflezl2rOa/fiU/VlDez5nHz6lPZvZ\nY4fjW7aEdu3CbN3yGVs2rmI3+wZrtpEDuuxF34G92GuvluTn+0c6+fk+MWnd2v/bqpU/frfd/NKy\npZ+HJi8P6pm7L+nUGFRERDJW0ElGQ6aXN4O99vLL52/exENuZb3jnfTrVMR9i9/hiy/8RHRr18Lq\n1Zv58x/u4fCvtpJPO9bRhnW045M385m+fC2t8vdnw4YQW7Y0LvZmzXyvm2bNti3RyYcZXHYZ/Pa3\ncRdPWlKiISIiGaG2S/GVVVUNGoW0ISO45m9dT6dOjs6dt+01vvQapn85ncHUMY/N1hDLfuTnsdmy\nBTZsgE2b/LJxo//3v//1j242b/Y/19T4ZcsW/9hn69ZtS+xDhWOOaXIxpR0lGiIikhGip5evVTu9\nvKtjFNJ4R3BdMmcOE3Yyj82U2bNh2jRatIA99/SL1C9NnhiJiIjsXDwT2R0/ZAjz6mkcUdcIro2Z\nx0YaRomGiIikvXgTgDGTJjGlqIi5odB3M8Y6YG6k++1VMd1vNbts4inREBGRtBdvAlA7guuyUaMY\nWFDAsA4dGFhQwLJRo+odJr6xtSCyc2qjISIiGSHeiewaO4KrZpdNLNVoiIhIRmjsY5C6NOSRRzy1\nIFI/DdglIiIZI9mjkIJGBm0qPToREZGMkYqJ7LIpyUgFPToREZGMpAQgMyjREBERkcAo0RAREZHA\nKNEQERGRwCjREBERkcAo0RAREZHAKNEQERGRwCjREBERkcAo0RAREQlAJo68HQQlGiIiIglSXV3N\n+NJS+hcWckbHjvQvLGR8aSnV1dWpDi1lNAS5iIhIAlRXVzO8pIQrq6qYEDXr67zp0xm+cGHOTsim\nGg0REZEEuGXsWK6MmloewIDB4TCjq6qYPG5cKsNLGSUaIiIiCbBkzhwGhcN1bhscDrNk9uwkR5Qe\nlGiIiIg0kXOO1jU11DfNmwGtampysoGoEg0REZEmMjM25uVRXxrhgI15eTk546wSDRERkQQ4fsgQ\n5oXq/lh9OhTihKFDkxxRelCiISIikgBjJk1iSlERc0Oh72o2HDA3FGJqURFXTZyYyvBSRomGiIhI\nAuTn5zOzvJxlo0YxsKCAYR06MLCggGWjRuVs11bQOBoiIiIJk5+fz4Rp02DaNJxzOdkmI5ZqNERE\nRAKgJMNToiEiIiKBUaIhIiIigVGiISIiIoFRoiEiIiKBUaIhIiIigVGiISIiIoFRoiEiIiKBUaIh\nIiIigVGiISIiIoFRoiEiIiKBUaKRBcrKylIdQlpQOWyjsvBUDtuoLDyVQ/IFkmiY2R5mdquZrTSz\nTWb2gpkdvZP9+5pZOGbZamb7BBFfttEfjqdy2EZl4akctlFZeCqH5Atq9tZ7gG7ABcBq4MfAM2ZW\n5JxbXc8xDjgYqP5uhXOfBxSfiIiIJEHCazTMbDfgLOBq59wS59wHzrnrgX8Dl+3i8LXOuc9rl0TH\nJiIiIskVxKOT5kAz4L8x6/8DnLCT4wx4zcw+NbN/mtlxAcQmIiIiSZTwRyfOuQ1mVg78r5ktB9YA\n5wMlwL/qOWw1cCnwCtAS+DmwyMx6Oedeq2P/3QCqqqoSHX5GWrduHZWVlakOI+VUDtuoLDyVwzYq\nC0/lsN1n527JuJ455xJ/UrNC4F6gL/AtUAm8BxQ75w5r4DkWAR86535Sx7bzgYcSFrCIiEjuucA5\n93DQFwmkMahzbgVwspntDnzPObfGzB4BVjTiNC8Bx9ezbR6+oelKYHNTYhUREckxuwEF+M/SwAVS\no7HDRczaAR8AY5xz9zTwmH8C651zZwcanIiIiAQmkBoNMxuIb9z5LtAVuBmoAu6PbP8d0KH2sYiZ\nXYGv7Xgbn2n9HDgZGBBEfCIiIpIcQY2j0Qa4EegAfAU8Boxzzm2NbN8P6Bi1fwtgMrA/sAl4A+jn\nnFscUHwiIiKSBEl5dCIiIiK5SXOdiIiISGCUaIiIiEhgUpZomFkfM5ttZqsik6gNjdne2sxuN7OP\nIxOzvW1ml8bs097M/mJmq81sg5lVmNlZMfu0M7OHzGydmX1tZnebWetkvMeGakBZ7GNm90e2bzSz\np8zsoJh9WprZdDP7wsyqzeyx2EnpzKyjmT0ZOcdnZnazmaVNstnUcoj8X99mZssj2z80s2lm9r2Y\n86R1OUBifidi9p9bz3nSuiwSVQ5mVmJmCyL3iXVmtsjMWkZtT+v7RILuEdlyv/y1mb1kZuvNbI2Z\n/d3MDo7ZJyH3QzM7KVJOm83sPTPbYVynVElEOZhZDzN72Mw+sm2fs6V1XKtJ5ZDKG0pr4DVgJH5C\ntVhTgYH4UUUPBW4Fbjez06P2+Qu+V8vpwGHA48CjZnZE1D4PA0VAP+AHwInAnQl9J023q7KYhe/z\nPAQ4EvgIP0nd7lH73Ip/f8Px73F/YGbtxsgf0FP4BsC9gZ8A/wPckNB30jRNLYf98Q2Nr8T/PvwE\nGAzcXXuCDCkHSMzvBABmNhrYGnueDCmLJpeDmZUAc4GngaMjy+1AOOo86X6fSMTvQ7bcL/sAfwSO\nBfoDecA/E30/NLMC4B/AAuAIYBpwt5mlS2/IeMvh8ajtPYHP8eNSdQMmATea2eW1OySkHJxzKV/w\nf/BDY9a9CYyNWfcKcEPU62r8yGbR+3wB/DTyc1Hk3EdFbR+EH61031S/74aUBf7GEAYOjVpn+KHd\na9/n9/Bzy5wZtc8hkeN6RV6fCtQAe0ftcynwNdA81e87EeVQz3nOxs+zE8rEcmhqWURuDB8C+9Rx\nnowqi3jLASgHJuzkvIdm0n2iCeWQdffLSIx7R+I+IfI6IfdD4CbgjZhrlQFPpfo9J6oc6jnP7cAz\nUa+bXA5pU0VahxeBoWa2P4CZnYz/g4oeyWwJ8MNIdZ+Z2Xn4uVIWRbb3Br52zr0adcwz+G8ExwYc\nf6K0xMf73SR1zv9P/5dtk9Qdjc/MF0Tt8y7+W01JZFVv4E3n3BdR556H74rcPajgE6gh5VCXtviB\n32q/vWZ6OUADyyLyzeZhYKSrezbkTC+LXZaDmX0f/7f+hZktiVSRLzKz6FGHS8js+0RD/zay9X7Z\nFh/jV5HXPUnM/bA3/v0Ts08J6SmecqhLm6hzQALKIZ0TjV/gB/n6xMy24Ku5RjrnlkTt80P8GBxf\n4v+o/oTP3j6IbN8XXy30HefH8vgqsi0TLMf/YtxoZm3NrIWZXQscgH9MANAe2OKcWx9z7Bq2vc99\nI69jt0NmlEVDymE7ZrY3MI7tq34zvRyg4WUxFXjBOfePes6T6WXRkHLoEvl3PP73YBB+7qUFZnZg\nZFum3yca+vuQdfdLMzP844EXnHPvRFbvS2Luh/Xt8z2Lat+TDppQDrHnOQ44l4bdMxtcDumcaJTi\ns+jTgWLgKuAOMzslap+J+OzrFHz2NgX4m5nt6tuYUfdzzrTjnPsWOAs4GP8HvwE/Wd1T+OfuO9PQ\n95n2ZdHYcjCzfOBJ4C3g+oZeJiHBBqwhZWG+seApwOh4L9P0SIPVwN+J2nvcn51zDzrnXnfOXYkf\ntfinu7hERtwnGvG3kY33yzvwbQtGNGDfRNwPrQH7pEKTy8HMDgOewD9mXLDDUTueg7rOU5egRgZt\nEjPbDd8oZZhz7unI6rfM7ChgDLDQzLrgG0Z1c84tj+zzppmdGFl/OfAZ/tl09LmbAe3YMUNLW5Gq\nzOLIh2cL59yXZrYUeDmyy2dACzP7Xkz2ug/b3udnwDExp24f+TcjyqIB5QCAme2Br9r7BjjLbRuR\nFrKgHKBBZXEy/tv8Ov9l5zuPm9li59wpZEFZNKAcVkf+rYo5tAroFPk54+8TuyqHbLxfmtntwGlA\nH+fcp1Gbmno//Czq3/Yx++yDfxS7panxJ0oTy6H2HN3wj0f+7Jy7MeYSTS6HdK3RyIsssdnSVrbF\n3CqyfWf7lANtIwlKrX74bGxZIgNOBudcdeQG0hXfLuOJyKYKfIOtfrX7Rro5dcK3dQFfFodHHifU\nGgisA94hg+ykHGprMv6JbwA6tI4/hKwpB9hpWdwI9MA3Bq1dAK4ALor8nDVlUV85OOdWAp/iG8FF\nOxjfSBay6D6xk9+HrLpfRj5chwEnO+c+itnc1PthVdQ+/djewMj6tNCEciiPWtcdWAjc55y7ro7L\nNL0cgmoBu6sF313rCHxXrDDwy8jrjpHtz+LnPOmL77b1P/h5UC6JbG8OvIdvyHQM/tvbVZGCHRR1\nnafwvVWOwU87/y7wl1S97zjL4uxIORTif6lWAI/GnOOOyPqT8NWiS4Dno7aHgNfx3fx64J9VrwF+\nm+r3n6hyAPYAluK7ARbis/DapbbXSdqXQ6J+J+o4Z2xvhbQviwT9bVyB700wHDgQ+C2wESiM2iet\n7xMJ+NvIpvvlHZH/zz4xf+O7xezTpPsh/nNnA77XxSH4Wp8tQP9Ul0ECy6E7vl3OgzHniO6N0+Ry\nSGUh9Y38wWyNWe6NbN8HuAf4OHJTeAe4IuYcBwJ/w1ePVgOvAufH7NMWmIHPVL8G7gJapfqXpJFl\n8Qt8Y6/NkV+aCcR0P8S3Hv8jvrtadaRc9onZpyO+P/SGyB/VTUQ+gNNhaWo5RI6PPbb2fJ0ypRwS\n9TtRxzm3smM38rQui0SVA3ANvgajGngBKInZntb3iQTdI7LlfllXOWwFLozaJyH3w0i5V+BrSP8F\n/DjV7z+R5YBvJF3XOT5IZDloUjUREREJTLq20RAREZEsoERDREREAqNEQ0RERAKjRENEREQCo0RD\nREREAqNEQ0RERAKjRENEREQCo0RDREREAqNEQ0RERAKjRENEREQCo0RDREREAvP/AYS7nDGPMDGl\nAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1063fba10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.plot(x,t,'ro')\n",
    "plt.plot(testx,testt,'b')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bendy!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Appendix: solving linear systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our expression for the optimal value of $\\mathbf{w}$ is:\n",
    "\n",
    "$$ \\mathbf{w} = \\left(\\mathbf{X}^T\\mathbf{X}\\right)^{-1}\\mathbf{X}^T\\mathbf{t} $$\n",
    "\n",
    "Above, we have computed this by writing code to compute the right hand side. This involves performing a matrix inversion *which is almost always something to avoid* (as the matrix gets bigger, this becomes time consuming and numerically innacurate).\n",
    "\n",
    "Instead, we can go back a step to this equation:\n",
    "\n",
    "$$ \\mathbf{X}^T\\mathbf{X}\\mathbf{w} = \\mathbf{X}^T \\mathbf{t} $$\n",
    "\n",
    "This is a system of linear equations (in general: $\\mathbf{A}\\mathbf{z} = \\mathbf{B}$) and we can solve this directly for $\\mathbf{w}$ using `numpy.linalg.solve`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "XX = np.dot(X.T,X)\n",
    "Xt = np.dot(X.T,t)\n",
    "w = np.linalg.solve(XX,Xt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[  4.55597856e+02]\n",
      " [ -4.43160485e-01]\n",
      " [  1.10151552e-04]]\n"
     ]
    }
   ],
   "source": [
    "print w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this system, it doesn't make any difference, but you don't have to get much bigger to see a difference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
